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ABSTRACT 

In this paper we present the analytic solutions for two test problems involving two-fluid mix- 
tures of dust and gas in an astrophysical context. The solutions provide a means of bench- 
marking numerical codes designed to simulate the non-linear dynamics of dusty gas. The 
first problem, DUSTYBOX, consists of two interpenetrating homogeneous fluids moving with 
relative velocity difference. We provide exact solutions to the full non-linear problem for a 
range of drag formulations appropriate to astrophysical fluids (i.e., various prescriptions for 
Epstein and Stokes drag in different regimes). The second problem, DUSTYWAVE consists of 
the propagation of linear acoustic waves in a two-fluid gas-dust mixture. We provide the ana- 
lytic solution for the case when the two fluids are interacting via a linear drag term. Both test 
problems are simple to set up in any numerical code and can be run with periodic boundary 
conditions. The solutions we derive are completely general with respect to both the dust-to- 
gas ratio and the amplitude of the drag coefficient. A stability analysis of waves in a gas-dust 
system is also presented, showing that sound waves in an astrophysical dust-gas mixture are 
Unearly stable. 

Key words: hydrodynamics — methods: analytical — methods: numerical — waves — ISM: 
dust, extinction 



1 INTRODUCTION 

Dust - from sub-micron-sized grains to centimetre-sized pebbles 
- is involved in many astrophysical problems. In particular, it pro- 
vides provide the materials from which the solid cores required for 
the planet formation process are built (see e.g. [Chiang & Youdin] 
|2010[ l. Dust grains are also the main sources of the opacities in 
star-forming molecular clouds, thus controlling the thermodynam- 
ics. Furthermore, observations are mostly sensitive to the dust - 
rather than the gas - emission. With the advent of the Spitzer and 
Herschel space telescopes our observational knowledge of dust at 
different wavelengths in young stellar and planetary objects has im- 
proved substantially. Millimetre and sub-millimetre observations 
will similarly be vastly improved with the arrival of ALMA that 
wiU achieve a spatial resolution < 0.1" at millimetre wavelengths 
( |Tumer & Wootten,2007) . 

Consequently, numerical simulations of astrophysical dust- 
gas mixtures are essential to improve our understanding of the sys- 
tems we will be able to observe. A dust-gas mixture is usually 
treated using a continuous two-fluid description, and a large class of 
numerical solvers have been developed. In an astrophysical context, 
two types of methods are generally adopted: grid-based codes (e.g. 
Fromang & Papaloizou 2006 , Paardekooper & Mellema|2006[ |Jo-| 
hansen et al.|2007||Miniati|201()^ or particle-based Smoothed Parti- 



cle Hydrodynamics (SPH) codes (e.g. Monaghan|1997[ [Maddison] 



let al.|2003|[Barriere-Fouchet et al.|2od5^ 



However, even with a continuous description of the mixture, 
the equations remain too complicated to be solved analytically for 
most problems, which presents a major difficulty for benchmarking 
numerical codes. Currently, the only known analytic solution in use 
is the solution for two interpenetrating homogeneous flows, given, 
e.g., by ,Monaghan & Kocharyan ( 1995 ) and Miniati (20101 for a 
linear drag regime and extended to one particular non-linear regime 
by |Paardekooper & Mellema| ( |2006| l. Knowing the analytic solution 
even for this simple case allows a precise benchmark of the various 
drag prescriptions that are appropriate in different astrophysical en- 
vironments (e.g., |Baines et al.|1965[ >. On the other hand, no usable 
analytic solution exists for the propagation of waves in a dust-gas 
mixture in a regime relevant to astrophysics, despite the rich liter- 
ature on the topic in the many other areas where dust-gas mixtures 
are of interest (for example in aerosols, emulsions or even bubbly 
gases, c.f. |Marble|1970||Ahuja|1973l|Gumerov et al.|1988[|Temkin| 
|1998[ l. Such solutions are of great interest as 1) they constitute a 
demanding test for a code's accuracy since small perturbations are 
easily swamped by numerical noise and 2) have to be correctly sim- 
ulated as they often appear in physical simulations. In the absence 
of such a solution, astrophysical codes (e.g. [Youdin & lohansen] 
[20071 pmia51[20T0l [BaT & Stone|[20T0) have generally been val- 
idated against the linear growth rates for the streaming instability 
( [Youdin & Goodman||200 5 1. Such a test problem is by definition 
limited to checking the growth rate of a given mode rather than val- 
idating against a full analytic solution. Another approach has been 
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to study numerical solutions for dusty-gas shock tubes ( [Miura &] 
|Glass|1982{|Paardekooper & Me llema 2006), where approximate 
solutions can be derived ^Miura& GlassI 19821 ) but again no com- 
plete analytic solution exists. 

In this paper we present the full analytic solutions for two spe- 
cific problems concerning two-fluid gas and dust mixtures in as- 
trophysics. The first, DUSTYBOX (Sec. |2|, is an extension of the 
interpenetrating flow solutions discussed above to the main drag 
regimes relevant to astrophysical dusty gases (i.e., Epstein and 
Stokes drag at different Reynolds and Mach numbers). The second, 
DUSTYWAVE (Sec.[3]( is the solution for linear waves in a dust-gas 
mixture, assuming a linear drag regime. 

Our aim is that these solutions will be utilised as standard tests 
for benchmarking numerical codes designed to simulate dusty gas 
in astrophysics. While it is beyond the scope of this paper to bench- 
mark a particular code using the two tests, the solutions we have 
derived were developed precisely for this purpose (for a new two- 
fluid SPH code that we are developing) and will be used to do so in 
a subsequent paper. 



V(justjin 

^diist-qiiad 
^dust.power 
Vdust.third 
^dust,mixed 
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2 DUSTYBOX: TWO INTERPENETRATING FLUIDS 

The first test problem, DUSTYBOX, consists of two fluids with uni- 
form densities pg and pd given a constant initial differential ve- 
locity (Avo — Vg,o — Vd,o). We assume that the gas pressure P 
remains constant. This test is perhaps the simplest two-fluid prob- 
lem that can be set up, for example by setting up two uniform 
fluids in a periodic box with opposite initial velocities. Thus it is 
a straightforward test to perform in any numerical code. Similar 
tests have been considered by Monaghan & Kocharyan (1995) for 
a single grain with a linear drag coefficient and by Paardekooper, 
|& Mellema| ( |2006t for one particular non-linear drag regime. How- 
ever, the simplicity of this test means that it can be used to test the 
correct implementation in a numerical code of both linear and non- 
linear drag regimes relevant to astrophysics, for which we provide 
the full range of solutions. 



2.1 Equations of motion 



The simplified equations of motion are given by: 



dvg 
' dt 
dvd 



pd = Kf ( Av) (vg - Vd) 



(1) 

(2) 



where momentum is exchanged between the two phases via the 
drag term (K being an arbitrary drag coefficient) and the function 
/(Av) specifies any non-linear functional dependence of the drag 
term on the differential velocity (i.e. / = 1 in a linear drag regime). 
In formulating ([TJ-l|2j it has been assumed that the effect of the col- 
lisions between the dust particles are negligible (i.e., no dust pres- 
sure or viscosity); that the dust phase occupies a negligibly small 
fraction of the volume (i.e., zero volume fraction: the estimated vol- 
ume fraction is ~ 10^^^ in planet forming systems); that the gas 
is inviscid; the two phases are in thermal equilibrium and that the 
only way for the two phases to exchange momentum comes from 
the drag term (that is, additional terms due to carried mass. Basset 
and Saffman forces have been neglected). 



Figure 1. Examples of the analytic solutions for the decay of the dust ve- 
locity in the DUSTYBOX test, assuming a dust-gas mixture with pg = 1, 
Pd = 0.01, i)d,o = 1, Vgfi = and K = 1 for the linear, quadratic, 
power-law (with a = 0.4), third order expansion (with 03 = 0.5) and the 
mixed (with 02 = 5) drag regimes. 



2.2 Analytic solutions 

Defining the barycentric velocity according to 

^. ^ PgVg.O + PdVd,0 

Pg + Pd 

and adding ([TJ and ([2j shows that the solutions to this equation set 
are of the form: 



Vg(i) = v* + 
Vd(i) = ' 



+ Pd 



Av(i), 
-Av(t). 



(4) 
(5) 



Pg + Pd 

The evolution of the differential velocity Av {t) depends on 
the drag regime. If the initial velocities of the two fluids have the 
same direction (say x), Eqs. |4](-l|5j reduce to two coupled scalar 
equations: 

Pd 



Vd,x (t) = 



K H V — (*) ■ 

Pg + Pd 

Avx (t) , 



+ Pd 



(6) 
(7) 

(8) 



where Avx is given by the differential equation 

dt VPg Pd/ 

The analytic expression for Av, (t) in five drag regimes 
f{Av) relevant to astrophysics in this particular configuration are 
given in Table [T] The linear solution (top row) holds for Epstein 
drag at low Mach number and Stokes drag at low Reynolds num- 
ber. A quadratic relation (second row) is relevant for Epstein drag 
at high Mach number and Stokes drag at large Reynolds numbers. 
Power-law drag occurs for Stokes drag at intermediate Reynolds 
numbers (in which case the exponent is given by a = 0.4). The 
third order expansion has been proposed for Epstein drag at in- 
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Drag type / Av^ {t) 

-kI 

Linear 1 oe 



Quadratic 



Pg Pd 



Power-law 



Tliird-order expansion 1 + 03 Ad^, 03 > 



Mixed 



V^rTa2AuJ, a2 > 



At)i,,oe 



l + asAj;^ ^ I 1-e 



-2K -!- + - 



sinh 




k Pg Pd y 




+ ^1 + a2 AiJ^ Q cosh 


[■< 


f 1. + 




cosh 


{" 


{- + - 

V Pg Pd y 




+ ^1 + a2AD2_osinh 




( ±. + 

Kph Pd) 





- 1 



Table 1. Expressions of Av (t) for several drag regimes / in a solid-gas mixture where the two phases have initially two different velocities. The pressure and 
densities of the medium are constant and the volume of the dust particles is neglected, e = +1 if Aiij;^o > and e = — 1 if Av^fl < 0. 



termediate Mach numbers ( [Baines et al.||1963] l. The mixed drag 
regime (bottom row) connects the linear and quadratic regimes for 
Epstein drag, used recently by Paardekooper & Mellema| (j2006 1. 

Finally, it should be noted that the stabiUty of the DUSTYBOX 
problem, though likely, has only been verified numerically. Proving 
stability with full generality is a difficult problem due to the non- 
zero mean velocities for each fluid — producing a dispersion re- 
lation that is a quadratic equation with complex coefficients. How- 
ever, it can be shown that the solution is stable for particular choices 
of K, vo, Cs and pQ. 



2.3 DUSTYBOX example 

As an example, the standard linear Epstein drag regime ( |Baines| 



et al.|196"5^ would correspond to K = pgCs / {pints) (where Cs is 
the sound speed, pi„t is the intrinsic density of the dust grains and 
s is the grain size) and / = 1. Thus, using the solution from Table 
[T] we would obtain the complete expression for the velocity in each 
phase according to 



+ Pd 



-Av^.oe VPS "dj , 



(9) 
(10) 



Examples of the solutions for the decay of the dust velocity in the 5 
different drag regimes for a typical astrophysical dust-gas mixture 
(i.e., 1% dust-to-gas ratio) characterised by pg = 1, pd = 0.01, 
"Ud.o = 1, ^g,o ~ and assuming K — 1 are shown in Fig.^ 
It may be observed, for example, that for this particular choice of 
parameters the quadratic and power law drag regimes (which would 
correspond to using a Stokes instead of an Epstein drag prescription 
in an accretion disc calculation) give less efficient relaxation of the 
dust phase to the barycentric velocity. 



3 DUSTYWAVE: SOUND WAVES IN A DUST-GAS 
MIXTURE 

The second test problem, DUSTYWAVE, consists of linear sound 
waves propagating in a uniform density two-fluid (dust-gas) 
medium with a linear drag term. Similar to the first problem, this 
test can easily be performed in 1, 2 or 3 dimensions in any numer- 
ical code using periodic boundary conditions. As previously, we 
provide analytic solutions for an arbitrary linear drag coefficient 
and/or dust-to-gas ratio. 

The setup consists of a sound wave propagating in the x- 
direction. We introduce the coefficient Cs, so that a small pertur- 
bation in the gas density 5pg is related to a small perturbation in 
the gas pressure 5P by the relation SP — c^Spg. Cb is thus the 
sound speed of the gas phase if no dust were present. 



3.1 Equations of motion 

For this system, the equations of motion are: 



Pg 



pd 



dt 



'■ dx 



dvd ^ ^ dvd 
dt dx 



dpg , dpgVg 



dt 



+ 



dx 



dpd ^ dpdVd 
dt dx 



-K(vg 



= +K{vg 



= 0, 



0. 



Vd) 



Vd) 



dP 

dx ' 



(11) 



(12) 



(13) 



(14) 



The assumptions made in obtaining these equations are identical to 
those discussed in Sec.|2] except for the additional term due to the 
gas pressure gradient. 
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3.2 Linear expansion 

We assume that the equilibrium velocities and densities of the fluid 
mixture are given by: Vg = Vd ~ 0, Pg = Pg,o and pd = pd.o- We 
then consider small perturbations and perform an acoustic linear 
expansion of Eqs. ljTTJ-l|T4]l: 



Pd.O 



dt 



+ 1 



9ug 

dvd 
dt 

dvg 
dx 



-K (ng — Vd) — c. 



= +K{vg~'Vd)., 



= 0, 



(Mpg 
dx 



(15) 



(16) 



(17) 



c»(5pd dvd _ 

-^^P"'"^ - °- ^^^^ 

As this system is linear, we search for solutions under the form of 
monochromatic plane waves. The total solution is a linear combi- 
nation of those monochromatic plane waves whose coefficients are 
fixed by the initial conditions. The perturbations are assumed to be 
of the general form 



Vi[kx- 



ot) 



Vd 

Spd 



i/gt , 

Dde^''''-'''\ 



(19) 
(20) 
(21) 
(22) 



where in general the perturbation amplitudes Vg, Vd, -Dg and Dd 
are complex quantities. 



3.3 Linear solutions 

Using l|19[l-l|22[) in l|T5j-l[T8j, we find that the resulting system ad- 
mits non-trivial solutions provided the following condition holds: 



1 

ILJ -\ 


1 


^ikc^ 
Ps,o 





1 

^td 










ikpgfi 





—iuj 








ikpd,o 





—iuj 



(23) 



Pg,o 



Pd.o 



where we have set tr, = — and td = ■ This condition pro- 

A K 
vides the dispersion relation of the system. 



3 I • 2 
UJ + lid 



1 



+ 



, 2 2 



. k Cg 

I 

td 



0. 



(24) 



This cubic equations admits three complex roots aj„=i,2,3 whose 
imaginary parts are always negative, ensuring the linear stability of 
the system (see Appendix[A|. This implies that the full solution of 
the problem consists of a linear combination of three independent 
modes that will take the form of exponentially decaying monochro- 
matic waves. For example, for the the gas velocity the solution will 
be given by 



Vg{x,t) 



+ 
+ 



g"iit [Vig J. cos {kx — LUiit) — Vig,i sin {kx — coirt)] 

g"2it [l/jg,! cos {kx — L02it) — V2g,i sin {kx — L02rt)] 



[Vs. 



g,r cos {kx — ujsit) — Vsg,! sin {kx — UJ-^rt)] ■ 

(25) 



where the subscripts r and i refer to the real and imaginary parts 
of the complex variables cji,2,3 and 14,1,2,3. The solutions for Vd, 
Spg and 5pd are of the same form, with the amplitudes replaced by 
the real and imaginary parts of Vd, Dg and Dd, respectively. 



3.4 Solving for the coefficients 

Obtaining the full analytic solution thus requires two steps: 

(i) Solving the cubic equation \24\ to determine the complex 
variables uui, 102 and lo^ for the 3 modes; and 

(ii) Solving for the 24 coefficients determining the amplitudes 
Vg, Vd, Dg and Dd for each of the 3 modes. 

Step i) can be achieved straightforwardly using the known ana- 
lytic solutions for a cubic equation (given for completeness in Ap- 
pendix |BJ. Since in general such solutions require a cubic equa- 
tion with real coefficients, it is convenient to solve for the variable 
UJ = —iy, which reduces Eq. \24\ to the form 



y -y 



1 1 

- + 



td 



,,22 

+ k c^y^ 



k Cg 

~tr 



= 0, 



(26) 



which has purely real coefficients, as required. 

Step ii) is less straightforward and consists of two substeps. 
The first substep is to constrain the amplitude coefficients using the 
8 constraints given by the initial conditions (i.e., the phase and am- 
plitude of the initial mode in the numerical simulation, which con- 
strain both the real and imaginary parts of the initial amplitudes). 
Although the solution can in principle be found for any given com- 
bination of initial perturbations to v and p for the two phases, the 
solutions we provide assume initial conditions of the form 



Vg{x,0) = Ug,o sin(fc3;), 
Vd{x,0) = Ud,o sin(fca;), 
Pg(a;,0) = pg,o +(5pg,osin(fca;), 
pd(a;,0) = pd,o + (5pd,o sin(fca::), 

giving the 8 constraints 

Vlg.i + V2g,r + V3g,r = 0, 

Vlg,i + V2g,i + V3g,i = -Vg,0, 

Vld.r + V2d,r + V3d,r = 0, 

Vld,i + 1^2d,i + V3d,i = -Vd,0, 

Dlg,r + D2g,r + D3g,r = 0, 

-Dlg,i + D2g,i + D3g,i = ~(5pg,0, 

Dld,r + D2d,r + D3d,r = 0, 

Did,i + D2d,i + Dad,! = — 5pd,o- 



(27) 
(28) 
(29) 
(30) 

(31) 
(32) 
(33) 
(34) 
(35) 
(36) 
(37) 
(38) 



The second substep is to determine the remaining 16 coef- 
ficients by substituting each of the expressions for the perturba- 
tions (Eq. |25]and the equivalents for Vd, 5pg and 5pd) and the 8 
constraints (|31^-(|38^ into the evolution equations (|15^-(|18^. The 
remaining analysis is straightforward but laborious, hence we per- 
form this step using the computer algebra system MAPLE. The re- 
sulting expressions for the 24 coefficients may be easily obtained 
in this manner, however their expressions are too lengthy to be use- 
fully transcribed in this paper Instead, we provide, for practical 
use, both the MAPLE worksheet and a Fortran (90) routing that 



^ Available as supplementary files accompanying the arXiv.org version of 
this paper. 
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Figure 2. Example analytic solution of the DUSTYWAVE test showing the 
propagation of a sound wave in a periodic domain in a two-fluid gas-dust 
mixture. The panels shows the time evolution (top to bottom, time as indi- 
cated) of the velocity in the gas (solid/black) and dust (dashed/red) respec- 
tively assuming Pg = Pd = 1 a gas-to-dust ratio of unity) and a drag 
coefficient K = I giving a characteristic stopping time of = 1/2. The 
solution with K = I shows efficient damping of the initial perturbation in 
both fluids. 

evaluates the analytic expressions for the coefficients (produced via 
an automated translation of the MAPLE output). The Fortran routine 
has been used to compute the example solutions shown in Figures|2] 
and [3] Note that, although the initial conditions are constrained to 
be of the form l|27|l-l|30^, the solutions provided are completely 
general with respect to both the amplitude of the drag coefficient 
and the dust-to-gas ratio. The examples we show employ a dust- 
to-gas ratio and drag coefficients that are typically relevant during 
the planet formation process. For this test it should be kept in mind 
that the solutions assume linearity of the wave amplitudes and do 
not therefore predict possible non-linear evolution of the system 
— for example the potential for mode splitting/merging or self- 
modulation, effects that are known to occur in multi-fluid systems. 

3.5 DUSTYWAVE examples 

Figure |2] shows a typical time evolution of the gas and dust ve- 
locities (solid/black and red/dashed lines, respectively) assuming a 
dust-to-gas ratio of unity (as would occur during the late stages of 
the planet formation process) and a drag coefficient K — 1, giving 
a characteristic stopping time of tj, — 1/[K{1/ pg + 1/ pd)] = 1/2. 
At t = 1 (bottom panel) it may be observed that the differential ve- 
locity between the two phases has been efficiently damped by the 
mutual drag between the two fluids. 

Figure|3]demonstrates how the characteristics of the solution 
change as the drag coefficient varies, again assuming a gas-to-dust 
ratio of unity, with the solution shown at t = 5. At low drag (top 
panel, K = 0.01, e.g. for dust in the interstellar medium) both the 
dust and gas evolve essentially independently over the timescale 
shown. Thus, the solution in the gas is simply a travelling wave 
with a sound speed close to the gas sound speed, while the dust 
retains its initial velocity profile. As the drag coefficient increases 
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Figure 3. Further examples of analytic solutions to the DUSTYWAVE test, as 
in Fig.|2]but showing the solution at f = 5 for a range of drag coefficients 
K = [0.01, 0.1.1.10, 100] (top to bottom, as indicated in the legend). At 
low K the waves are essentially decoupled in the two fluids (top panel), 
while an intermediate drag coefficient produces the most efficient damping 
(middle panels). At large K (bottom panels), although the differential ve- 
locities are quickly damped the overall amplitudes decrease more slowly 
since the waves tend to move together 



to unity (top three panels), the solution tends towards the efficient 
coupling that occurs ai K ^ 1 (for this gas-to-dust ratio) shown 
in Fig. [2] which represents the "critical damping" solution where 
both the gas and dust velocities relax to zero. As the damping is in- 
creased further (bottom two panels) the damping of the differential 
velocities occurs in a fraction of a period, implying that, although 
the differential velocity between the fluids is quickly damped, the 
removal of kinetic energy is less efficient since the two waves es- 
sentially evolve together, relaxing slowly — but in tandem — to 
zero. 



4 CONCLUSION 

In this paper we have provided the analytic solutions to two prob- 
lems involving two-fluid astrophysical dust-gas mixtures in order 
to supply a practical means of benchmarking numerical simula- 
tions of dusty gas dynamics. The test problems are simple to setup 
for both particle and grid-based codes and can be performed us- 
ing periodic boundary conditions. A summary of the setup for each 
problem is given in Table|2] It may be noted that both solutions are 
completely general with respect to both the amplitude of the drag 
coefficient and the dust-to-gas ratio. The subroutines for comput- 
ing the DUSTYWAVE solution are provided as supplementary files 
to the version of this paper posted on arXiv.org. The DUSTYWAVE 
solution has also been incorporated into the SPLASl|^visualisation 
tool for SPH simulations ^ffice]2007j. 



^ httpV/users.monash.edu.au/^dprice/splash/ 
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Test problem Initial conditions 



Boundary conds. 



DUSTYBOX Vg = 


vor 






Periodic 


Vd = 


~vor 








Pg = 


Pg,0 








Pd = 


Pd,0 








Cs,0 = 


= const 








DUSTYWAVE Vg = 


t^g^O sin(k ■ r)f 




Periodic 


Vd = 


fd,o sin{k ■ r)f 






Pg = 


Pg,0 + 


Spg o sin(k ■ 


r) 




Pd = 


Pd,o + 


(5pd,o sin(k 


■r) 




Cs,0 = 


= const 









Table 2. Summary of the setup for each of the test problems. 



(i.e. = uji > 0), the solution diverges with time and the system is 
unstable. 

Physically, the evolution of the solid-gas mixture given by 
Eqs. l|15[l-l|18[l is characterised by three intrinsic time scales: the 
time tp = (fccs)~^ required for the gas pressure to equilibrate the 
gas phase, the time (defined according to t^^ = + t^^) 
required by the drag to relax the centre of mass of the fluid (see 
Sec.[2]( and td, the time required for the dust to force the gas evolu- 
tion. Thus, two ratios of these independent timescales are sufficient 
to fully describe the evolution of the system. Defining the following 
dimensionless quantities: 



Y = ytp, 



(Al) 
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td 



td' 



(A3) 



giving Eq. l|26[) in the form 



Po {Y) ^Y^ -Yh + Y -ed^O. 



(A4) 
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APPENDIX A: STABILITY OF LINEAR WAVES IN A 
DUST-GAS MIXTURE 

The roots of Eq. \26\ can be three real single roots or one real roots 
and two conjugated complex roots (plus all the degenerated cases). 
Noting y = + iyi, we see e.g. from Eq. \25) that if j/r < 



Determining the sign of the real part of the roots of Eq. \AA) 
specifies whether the system is unstable or not. For this purpose, 
we introduce: 



P = 



e-" 1 
27 + 2 ' 



D 



' 27 



_e eed 

108 6 



1 

27' 



(A5) 



(A6) 



(A7) 



The roots of Eq. \M) are denoted r_i, ro and r+i. We have the 
following three cases: 

• Case(i): P < 0, 

• Case (ii): P > and D < 0, 

• Case (iii): P > and D > 0, 

In cases (i) and (ii) one root is real (denoted ^) while the two 
remaining roots are the complex conjugates of each other (denoted 
a ± iP). Factorising Eq.|A4|using the three roots gives the relations 



p + 2q = e, 



2fia + (a^ + /3^) 



= 1, 



Ed- 



(A8) 



(A9) 



(AlO) 



The last equation implies fj, > 0. Combining i |A8[ >- (aTo1 > gives an 
equation for a of the form 



r / \ 3 2 a / 2 , ,\ — Ed) 

/(a) = a -ea + ^ (e + 1) ^ 



0, 



(All) 



which admits only one or three positive roots provided e — td ~ 
^ > 0. Indeed, this is the case since we have / (0) < 0, /' (0) > 

0, lim / (a) and a positive X axis value Qc = ^ for the centre of 

a — ^ + oo o 

symmetry of the cubic function. Therefore, in cases (i) and (ii), the 
real part of the complex roots is positive and the system is stable. 
In case (iii), the three roots are real. To determine their signs. 



Analytic solutions for two-fluid dust- gas problems 7 



we calculate the Sturm polynomials of Po {Y): 



Po{Y) = 

Pi(Y) = 

P2{Y) = 

PsiY) = 



3Y^ - 2eY + 1, 



2 2e" 

3 ^ ~9" 



243D 

- 3)' 



(A12) 
(A13) 

(A14) 
(A15) 



As the three roots are real, D > and P3 > 0. We can then apply 
the Sturm theorem to determine the number of positive roots. Us- 
ing V (Y) to denote the number of consecutive sign changes in the 
sequence [Pq (Y) , Pi (Y) , P2 {Y) , P3 (Y)], the number of posi- 
tive roots of Po (Y) is given by V (0) — lim V (Y). Thus, if 

ed — I < 0, the three roots are positive. However, if — f > 0, 
only one root is positive and the two remaining ones are nega- 
tive. However, ea has to be smaller than the larger positive root 
of D(ed) = for D in order to remain positive (and thus, for the 
system to be unstable), which is never the case if 9ed > e > \/3. 

Such a solid-gas mixture would thus be unconditionally sta- 
ble. However, if e is large enough (which means physically, that 
the damping due to the drag occurs faster than the equilibrium pro- 
duced by the pressure) and the ratio — is sufficiently large (i.e., a 

e 

large dust-to-gas ratio and thus, an efficient forcing of the gas mo- 
tion from the dust), and instability may develop when additional 
physical processes are involved (an example being the streaming 
instability that occurs in a differentially rotating flow, c.f. |Youdin| 
|& Goodnian|2005) . 



where 

1 • , / Q \ 

t — - arcsmn — ^^=^= 

Case ii) P > and D < 0, 



—a u + V /—u — V 



where: 

u — 



Case iii) P > and P» > 0, 

/ 2nn + arccos ( -3= 
r„ = ^ +2x/P,os 




(B7) 

(B8) 
(B9) 

(BIO) 
(BID 

(B12) 



withn = 0,±1. 



APPENDIX B: REAL AND IMAGINARY PARTS OF THE 
ROOTS OF A CUBIC WITH REAL COEFFICIENTS 

The cubic equation given by Eq.|26]can be solved using the known 
analytic solution to a cubic equation, though for this problem we 
require both the real and imaginary components of all three solu- 
tions. We consider the following normalised cubic equation with 
respect to the variable x: 

f (x) ^ x'^ + ax"^ + bx + c, (Bl) 

where a, b and c are real coefficients. We introduce the quantities 
P, Q and D, given by: 

(B2) 

^ 6 2 27' 
and 

D^P^-Q'^. (B4) 



Denoting the roots of Eq.l |Bl^ by r_i, ro and r+i, the solutions 
can be obtained by considering the following three cases: 

Case i) P < 0, 

ro = ^ + 2^/^ sinh (t) , (B5) 
r±i = ^ - sinh (t) ± iV-SP cosh (t) . (B6) 



